import numpy as np
import statsmodels.api as sm
import pandas as pd

# 生成数据
X = np.random.rand(100, 1) * 10  # 特征
y = 2.5 * X.squeeze() + np.random.normal(0, 2, 100)  # 目标变量 + 噪声

# 添加截距项并拟合模型
X = sm.add_constant(X)
model = sm.OLS(y, X).fit()

# 输出回归结果
print(model.summary())

# 残差诊断（可视化）
import matplotlib.pyplot as plt
residuals = model.resid
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].scatter(model.fittedvalues, residuals)
ax[0].axhline(y=0, color='r', linestyle='--')
ax[0].set_title("Residuals vs Fitted Values")
sm.qqplot(residuals, line='s', ax=ax[1])
plt.show()
